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INTRODUCTION 


I . 

Research  in  the  time  period  covered  by  this  report  has  been 
concentrated  primarily  on  the  subject  of  linear,  space-’var lent  restora- 
tion of  blurred,  photon  limited  images.  However,  some  effort  has  also 
been  devoted  to  extending  our  previous  work  on  photon  limitations  in 
coiTipensated  imaging,  particularly  an  extension  of  the  results  presented 
in  iXADC  report  mOC-TR-77-165. 

In  section  II  v^e  discuss  the  extensions  mentioned  above.  Section 
III  briefly  describes  some  unsuccessful  attempts  at  formulating  the  space- 
variant  filtering  problem  in  continuous  notation.  Section  IV  deals  with 
a more  successful  discrete  formulation  of  this  problem.  Finally,  section  V 
outlines  the  numerical  computations  now  being  undertaken. 

I I . FURTHER  PERFORMANCE  PREDICTIONS  FOR  A PRE-  AND 

POST-COMPENSATED  IMAGING  SYSTEM 

In  two  previous  technical  reports  (RADC-TR-76-382  and  RADC*TR-77"165) 
the  photon- li  mi  ted  performance  of  a sped  Me  compensated  imaging  system 
was  analyzed.  The  primary  results  of  this  analysis  were  two  figures 
showing  restored  bandwidth  aIT  and  quality  factor  0 ds  a function  cf  the 
total  number  N of  photoevents  intercepted  by  the  system  for  a variety 
of  combinations  of  ore-  and  post-processing.  Figure  8 of  RADC-TR-76-382 
showed  that  the  restored  bandwiuth  of  the  pre-  and  post -compensated 
system  reached  the  d i f f ract ton- 1 imi ted  bandwidth  when  approximately  10^ 
photoevents  were  intercepted  by  the  system,  this  conclusion  being  valid 
for  a point-source  object. 

In  the  derivation  of  these  curves,  certain  assumptions  about  the 
syste^m  were  made.  These  included  assumptions  that  the  atmospheric 
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coherence  diameter  r^  was  lOcm,  that  the  system  employed  a shearing 

Interferometer  with  317  subapertures,  and  that  the  ratio  of  Image 

4 

Integration  time  to  wavefront  sensor  integration  time  was  10  . In 
addition  It  was  assumed  that  the  system  employed  a fixed  splitting 
ratio  such  that  90^  of  the  incoming  light  was  sent  to  the  wavefront 
sensor  and  10'^  was  sent  to  the  Image  detector.  The  90^  splitting  ratio 
was  found  to  be  nearly  optimum,  although  the  maximum  is  very  broad. 

Concerns  wf»re  expressed  by  members  of  the  Review  Group  that  the 
number  of  photoevents  required  in  the  image  {\0%  x 10^  " 10^  photoevents) 
might  saturate  the  detector.  Examination  of  the  calculations  showed 
that  the  requirement  for  10^  total  photoevents  was  dontinated  by  the  flux 
required  for  satisfactory  operation  of  the  wavefront  sensor,  and  it  was 
speculated  that  the  use  of  s.-naller  splitting  ratios  (less  than  I0^^  to  the 
’mage)  would  reduce  the  dynamic  range  required  of  the  detector  without 
seriously  degrading  the  quality  of  the  final  restored  image. 

The  correctness  of  this  speculation  was  investigated  by  setting 
the  performance  of  the  wavefront  sensor  at  certain  fixed  levels,  and 
calculating  restored  bandwidth  and  intage  quality  as  a function  of  the  total 
number  of  photoevents  In  the  detected  image.  This  latter  quantity  is 
represented  by  the  symtx>l  N.. 

The  results  of  these  calculations  are  presented  in  Fig.  I.  Part 

(a)  shows  the  restored  bandwidth  aTT  (cycles/m  rad)  vs.  N.  , while  part 

(b)  shows  the  quality  factor  Q vs.  N..  In  both  cases,  three  curves  are 

presented,  one  for  each  of  three  levels  of  performance  of  the  wavefront 
sensor.  The  parameter  represents  the  residual  rms  wavefront  error 

fmeasurement  noise  plus  fitting  error)  and  takes  on  the  values  0.5  radians, 
1.5  radijns  and  4.5  radians.  These  numbers  correspond  to  0.08,  0.24  and 
0.72  waves  of  rms  error. 
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Figure  I:  Restored  bandwidth  (a)  and  quality  factor 

(b)  vs.  average  number  of  im.ige  photoevents 
for  fixed  levels  of  performance  of  the  wave 
front  sensor. 


As  can  be  seen  from  part  (a)  of  the  figure,  when  o « 0.5  radians, 
the  restored  bandwidth  is  essentially  indistinguishable  from  the  diffraction 
limited  bandwidth,  even  for  as  few  as  100  photoevents  in  the  image. 

For  0 • I.S  radians,  about  10  image  photoevents  are  needed  to  achieve 
the  diffraction-limited  restored  bandwidth.  Finally,  when  o • A. 5 
radians,  the  performance  of  the  wavefront  sensor  is  so  poor  that  enormous 
numbers  of  image  photoevents  (considerably  more  than  10^^)  are  needed  to 
achieve  diffraction-limited  performance. 

Hopefully,  this  set  of  calculations  will  provide  some  indication 
of  the  numbers  of  image  photoevents  that  will  be  required  to  achieve 
satisfactory  quality  in  the  final  restored  images. 

III.  CONTINUOUS  FORHUIATIOH  OF  THE  SPACE-VARI ANT 

RESTORATION  PROBIEH 

All  previous  analysis  of  the  performance  of  pre-  and  post -compensated 
imaging  systems  has  assumed  that  the  final  post-detect ion  restoration 
filter  is  a linear  space* invar iant  Wiener  filter.  This  type  of  filtering 
is  known  to  be  optimum  only  when  the  signal  and  noise  are  uncorrelated 
stationary  random  processes  and  when  the  noise  statistics  are  gauss ian. 

In  the  regime  of  photon- I imi ted  imaging,  the  noise  is  non-gaussian  and 
s igna I -dependent , and  it  is  therefore  reasonable  to  assume  that  iome  form 
of  linear,  space-variant  filtering  or  rionl  inear  filtering  will  perform 
better  than  the  simple  filtering  used  in  earlier  analyses. 

The  change  from  linear  space- invar iant  filtering  to  more  complex 
filtering  strategies  has  associated  with  it  a cost  in  computational  com* 
plexity  and  proces»ir>g  time.  It  is  natural  therefore  to  Inquire  as  to 
how  much  image  quality  is  gained  by  these  more  complex  methods,  and  at 
what  price. 
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In  this  report  we  consider  only  linear,  space-variant,  least-mean- 
square  filtering.  Such  filters  depend  on  certain  average  properties  of 
the  class  of  images  anticipated,  but  do  not  depend  on  the  particular  image 
detected  on  any  given  trial.  Techniques  which  perform  a filtering  oper- 
ation that  depends  on  the  particular  image  detected  are  generally  non- 
linear, and  will  be  considered  in  later  work  under  this  contract. 

Figure  2 illustrates  the  nature  of  the  lea$t-square-f i I ter ing 
problem  of  interest  to  us  here.  For  simplicity  a one-dimensional  space 
variable  x is  used  throughout.  The  object  radiance  distribution  o(x) 
is  assumed  to  be  a random  process.  In  general  the  statistics  o(x)  may 
be  non- s tat  ionary;  hence  its  autocorrelat ion  function 
a function  of  the  two  space  coordinates  and  X2  rather  than  just 
their  difference. 

The  object  radiance  distribution  is  subjected  to  a linear  blur, 
described  by  an  impulse  response  or  point- spread  function  that  is  in 
general  space-var iant , and  is  represented  by  b(xpX2)>  The  result  is  a 
classical  image  irradiance  i (x)  incident  on  the  detector,  where 

i(x)  • I b(x,c)o(c)(<C  (0 

Through  Che  detec c ton  process,  a detected  image  d(x)  is  generated. 

As  tinp)  ied  by  the  semi -c  lass  teal  theory  c>f  photodetection,  d(x)  is  a 
doubly  stochastic  Poisson  impulse  process,  with  space-variant  mean  X(x) 
related  to  the  classical  image  intensity  through 

Mx)  - 2I  i(«)  (2) 

hv 

where  ^ is  the  quantum  efficiency,  T the  integration  time,  h is 
Planck's  constant,  and  v is  the  mean  optical  frequency.  The  detected 


- 5 - 


Image  d(x)  is  passed  through  a linear,  space-variant  restoration  filter 
with  impulse  response  h(xj,X2)  » yielding  a restored  Image  r(x)  given 
by 

aa 

r(x)  - j h(x,n)d(n)dn  (3) 

•‘Od 

The  impulse  response  h(xpX2)  is  chosen  to  minimize  the  mean-square 
difference  between  the  restored  Image  r(x)  and  an  ideally  filtered 
object  o(x).  In  general,  the  ideal  filter  may  be  space-variant  with 
impulse  response  s(X|,X2)* 

At  this  point  it  is  natural  to  inquire  as  to  what  set  of  conditions 
can  lead  to  a I east -mean- square  filter  which  is  space-variant  rather  than 
space- invar iant . To  answer  this  question,  v/e  must  first  present  some 
theoretical  results  which  can  be  derived  without  much  difficulty,  and 
which  are  presented  here  without  proof.  The  impulse  response  h(x,x'^) 
of  the  restoration  filter  which  achieves  least  mean-square  error  is  the 
solution  of  the  Integral  equation. 

m 

h(x,x")R^(x'  ,x'^)dx"* 

- I s(x.x")R^^(»'.x")d*"  W 

Here  s(x,x")  is  again  the  known  impulse  response  of  the  ideal  filter, 
while  and  are  also  kno%^  ^unctions,  given  by 

R^(x‘,x")  - £id(x')d(x")] 

T 2-2 

- ».(*•. X")  (5) 

hv  h V ’ 
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RdoU-.x^') 


E[d(x')o(x")] 


- (x'.x")  (6) 

hv 

where  i(x')  Is  the  mean  value  of  the  Image  Irradlance  at  , 
averaged  over  the  ensemble  of  possible  objects,  R|(x^,x")  Is  the  auto- 
correlation function  of  the  Image  Irradlance,  and  R,  (x',x^^)  Is  the 

10 

cross-correlation  function  of  the  Image  Irradlance  and  the  object  radiance. 
More  specif  leal ly, 

•e 

I (x^)  ■ f b(x^,x^^)o(x^^)dx'‘^  . (7) 


R,(x-,x'')  - Ell(x')l(x'01  - II  b(x',Ob(x",n) 


Ro(E.n)<JC<in  . 


(8) 


Rjg(x",x"')  - E[1  (x')o(x'")l  • J b(x-,{)R^(5,x")d{  , 

(9) 

where  o(x'^)  Is  the  mean  object  radiance  at  x'^  » and  R^(c»n) 
autocorrelation  function  of  the  object  radiance  distribution. 

Examination  of  this  somewhat  bewildering  array  of  results  leads 
one  to  the  conclusion  that  the  Icast-mean-square-error  filter  will  be 
spacc“var iant  if  any  one  of  the  following  conditions  hold: 

(1)  The  mean  value  o(x)  of  the  object  radiance  distribution  Is  not 
constant  (l.<  . , is  indeed  a function  of  x)  , the  expectation 
being  over  the  entire  ensemble  of  possible  objects; 

(2)  The  autocorrelation  function  R^(C,n)  of  the  object  Is  non- 
statlonary  (i.e.,  depends  on  both  C and  n , rather  than  just 
their  d i f f erence) . 
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(3)  The  impulse  response  b(xj,x2)  space  variant  (i.e., 

non  i soplanat ic) ; 

(M  The  impulse  response  s(Xj,X2)  of  the  ideal  filter  is  space 

variant. 

We  now  examine  each  of  these  conditions  in  the  context  of  the 
compensated  imaging  problem.  Consider  first  the  mean  value  o(x) 
of  the  object  radiance  distribution  over  the  object  ensemble.  All  space- 
objects  are  of  course  spatially  bounded  (i.e.,  of  finite  extent), and 
for  this  reason  we  could  argue  that  o(x)  is  always  a function  of  x. 
However,  in  any  real  measurement  we  look  in  on  the  object  process  with  a 
finite  measurement  window,  and  it  is  possible  that  over  this  measurement 
window  the  expected  object  radiance  is  constant.  On  the  other  hand,  it  is 
not  difficult  to  identify  situations  in  which  the  mean  object  radiance 
would  be  a function  of  position  within  the  measurement  window.  An  example 
would  be  when  the  object  of  Interest  is  far  smaller  than  the  measurement 
window,  but  due  to  tracking  errors  its  position  within  the  window  obeys 
some  non-uniform  probability  density  function.  In  such  a case  o(x) 
is  definitely  a function  of  x , and  therefore  the  least-mean-square 
restoration  filter  must  of  necessity  be  space-variant. 

Turning  to  condition  (2),  we  would  also  expect  the  object  auto- 
correlation function  R^(C»n)  to  be  non-stii.  ionary  under  the  conditions 
described  above.  Alternatively,  even  if  there  are  no  tracking  errors,  and 
the  object  has  constant  mean  radiance  in  the  field  of  interest,  we  might 
expect  the  central  portion  of*  the  object  (i.e.,  the  main  body)  to  have 
different  textural  appearance  than  the  extremities  (appendages  consisting 
of  arms,  booms,  etc.),  and  hence  a non-stat ionary  autocorrelation  function 
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might  be  required  to  properly  describe  the  second-order  statistics  of 
the  object.  Again  a space-variant  restoration  filter  would  be  required 
for  optimum  performance. 

As  for  condition  (3)i  space  variance  of  the  blur  Is  anticipated 
when  the  size  of  the  object  exceeds  the  size  of  the  Isoplanatlc  region 
of  the  compensated  Imaging  system.  Under  such  conditions^  space-variant 
restoration  Is  again  required  If  the  smallest  possible  mean-square  error 
Is  to  be  achieved. 

Finally,  condition  (A),  space-variance  of  the  Ideal  filter,  can 
generally  be  Ignored.  In  most  cases  of  Interest  It  Is  appropriate  to  use 
as  an  ideal  filter  the  diffraction  limited  telescope  Itself.  Such  a filter 
Is  space  Invariant,  and  as  a consequence  It  Is  always  safe  In  these 
problems  to  assume  that  the  Impulse  response  s depends  only  on  the 
difference  of  Its  two  variables. 

At  this  point*  we  have  Introduced  arguments  that  demonstrate  that 
In  most  practical  applications  of  compensated  imaging,  the  Ideal  post- 
detection restoration  filter  will  be  space  variant.  It  is  quite  another 
matter  to  say  precisely  what  the  optimum  filter  Is,  and  to  calculate  its 
performance.  The  chief  difficulties  Involved  In  this  task  are: 

(1)  Specification  of  reasonable  models  for  o(x)  , R^(x,x')  and 

b(x,x^)  ; and 

(2)  Solution  of  the  Integra!  equation  {k)  for  these  models. 

Several  attempts  at  finding  a solution  for  various  specific  cases 
have  not  met  with  success.  Case$  examined  Included  an  object  which  is 
a 'Stfindowed"  stationary  process  (and  therefore  nonstat lonary) , and  an 
object  with  a nonstationary  autocorrelation  function  of  the  form 
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Since  even  these  idealistically  chosen  cases  proved  unsolvable,  it  was 
decided  that  the  continuous  approach  would  have  to  be  abandoned  in  favor 
of  a discrete  approach.  For  a discrete  formulation  of  the  problem,  the 
necessity  of  solving  an  integral  equation  to  find  the  required  impulse 
response  h is  replaced  by  the  necessity  of  solving  a matrix  equation. 

The  latter  is  far  easier  to  do,  in  general.  We  therefore  turn  to  a 
discussion  of  the  discrete  approach  to  the  problem  at  hand. 

1 V . DISCRETE  FORMULATION  OF  THE  SPACE-VARIANT 

RESTORATION  PROBLEM 

In  this  section  we  treat  the  same  problem  addressed  in  section  111, 
but  using  a discrete  representation  of  the  various  physical  quantities 
involved.  We  begin  with  a discussion  of  the  discrete  representation 
itself  and  the  physical  meaning  of  the  quantities  involved. 

(a)  Discrete  Notation 

The  continuous  object  radiance  distribution  o(x)  is  replaced  by 
a column  vector  o consisting  of  uniformly  spaced  samples  of  o(x)  , it 
being  assumed  that  the  object  is  approximately  bandlimited  and  can  there* 
fore  be  subjected  to  sampling  without  appreciable  loss  of  information.  In 
the  case  of  a two  dimensional  object  o(x,y)  , the  vector  o Is  con- 
structed by  lexicographic  order in?^ of  the  samples.  For  an  object 
consisting  of  ML  elements,  the  vector  o contains  M ■ KL  elements. 

The  blurring  operation  to  which  the  object  radiance  is  subjected 
can  be  represented  by  an  matrix  (B].  Such  a matrix  can  represent 

a space- invar i ant  or  a space-variant  blur,  depending  on  the  form  of  its 
elements.  The  elements  of  [B]  are  samples  of  the  impulse  response 
see  the  Appendix  for  a detailed  definition  of  this  term. 
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b(X|,X2).  The  classical  Intensity  of  the  image  falling  upon  the  photo- 
detector is  represented  by  a column  vector  given  by 

i ■ [B]o  (II) 

The  detector  is  assumed  to  consist  of  a discrete  array  of  elements, 

each  of  which  produces  a photocount  for  the  observation  interval  T used. 

Thus  the  detected  signal  is  represented  by  a column  vector  ^ of  length 

H ; it  is  assumed  that  the  n^^  element  d of  d is  a Poisson  random 

n - 

variable,  with  mean  given  by 

\ « i (12) 

n — n 

hv 

where  A is  the  area  of  a detector  element,  and  i is  the  n^^ 

n 

element  of  Equivalently,  we  write 

X . ■ ^[Blo  . (13) 

hv  hv 

where  the  n^^  element  of  X is  X . 

— n 

The  detected  signal  ^ is  subjected  to  linear  filtering  with  the 
aim  of  restoration.  The  restoration  filter  is  represented  by  a matrix 
(H] , the  elements  of  which  are  samples  of  the  (possibly  space-variant) 
impulse  response  h(x^X2)<  The  filter  (H]  is  to  be  chosen  to  minimize 
a certain  measure  of  error  between  the  restored  image  r_  given  by 

r - [H]d  (III) 

and  an  "ideally  filtered"  object  given  by 

5 ■ [S]o  (15) 

Here  [S]  is  an  matrix  of  samples  of  the  impulse  response 

s(X|,x2)  of  the  ideal  filter.  The  use  of  a " over  a matrix,  here  and 
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elsewhere  in  what  follows,  implies  that  this  matrix  has  been  normalized 
such  that 

I [S]^^  - I (16) 

m n 

where  [S]  ^ indicates  the  (n,m)  element  of  (S). 
nm 

We  define  an  error  column  vector  £ given  by 

tS]o  - [Hj^  (17) 

The  restoration  filter  [H]  will  be  chosen  to  minimize  the  mean-squared 
norm  of  the  vector  £.  Equivalently  we  choose  [H]  to  minimize  the 
quantity 

* " j ■ 

where  t signifies  a transpose  operation,  E is  an  expectation  operator, 
and  Tr(  ) signifies  the  trace  of  a matrix  (i.e.,  the  sum  of  the  elements 
along  the  main  diagonal). 

(b)  Minimizing  Hean-Squared-Error 

Substituting  Eq.(l7)  in  Eq.(l8),  the  mean-squared  error  to  be 
minimized  becomes 

*■  E}Tr(cc‘)} 

- E|Tr(([S]o-  (H]d)([Sjo-  (H]d)')l 

- E|Tr([H]d  d‘(H)‘  - 2(S)od‘(H]‘*  (S)o  o‘(S]*)|  (19) 

The  expectation  can  be  performed  in  two  steps.  First  we  assume  o is 
fixed  and  take  the  expectation  over  only  the  Poisson  statistics 

E^(«)  - E{*|  2.) 

Later  we  take  a second  expectation  over  the  statistics  of  o.  Thus 
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(20) 


EpjT'-(£e‘)|  ■ Tr([H]Epjd  d*j[H]‘ 

- 2[S]o  Ep|d*|[Hl*  + [S)o  o‘[S]*) 

Because  £ is  Poisson  distributed*  with  statistically  independent 
components  (assuming  o is  fixed)*  it  follows  that 


I t)  \ ® t 

EpddM  - X +ii 

b \ 

where  the  notation  is  defined  by 


0 \ 


'1  0 


0 **a 


It  Is  convenient  at  this  point  to  normalize  the  matrices  and  vectors 

for  future  expressions.  Again*  a matrix  [A]  is  normalized  to  produce 

M 

[A]  by  factoring  out  a constant  C - max  7 [A]  from  all  the  columns. 

m . nm 

n«l  H 

A vector  a is  normal ized  to  produce  a by  factoring  out  I a from 

M " n-I  " 

all  elements.  As  a consequence*  I i ■ 1.  This  norjaal tzaiion  assumes 

^ I ^ 

n»l 

that  the  amount  of  light  intercepted  by  and  the  net  energy  loss  of  the 
optical  system  does  not  change  as  a point  source  moves  about  the  object 


field. 


Now,  assuming  that  the  energy  loss  is  approximately  constant  over 


the  field  of  view*  (13)  can  be  written 

X - r l»lo 


- Ill  - 


where  we  have  used 


r 


A hv 

b nTA 
o 


b " y [B]  (assumed  identical  for  ail  m) . 
0 ^ nm 

n 


Using  (2^)  in  (22),  we  obtain 


E jd  d*j 

t>|-~  I 


^ (B)o  i+ 


In  addition 


E 'd^ 
P|- 


Finally,  substitution  of  (26)  and  (27)  in  (21)  yields  the  conditional 
expectation 

EpIfrUeS}  - Tr||Hl^  ^ ^(B]o  ^ + J_  (B]o  o‘(B)^yH]‘ 


i[iloo‘(B]‘[H]'  ♦ [i]o  o'[il‘ 


The  remaining  task  in  calculating  the  desired  error  measure  Is 
to  average  over  the  object  ensemble.  Applying  such  an  average  to  Eq.(28) 


yields 


6 - TrLllN  '^(B)5  ° + HBjMlHl' 

\ Uo  \j  ^ \ 

■ 2k  N^(S1(«  HB]*[H]‘  ♦ k^  N^[S)[*  l[S)*j 

O Of 
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where  the  expected  number  of  photoevents » N , in  the  image  is 


o)k  ^ ''n 
n 


(30) 


in  addition,  the  definitions 


o • E (o) 
— o 


£}/  I 


l«oJ 


E o o 


(31) 


have  been  used.  This  is  net  a totally  unusual  normalization  for  ] 

o 


since 


where 


I ((B]o) 

n 


- 1.2 


of  » variance 


ii  *") 


(32) 


(33) 


Thus,  if  the  variance  of  the  number  of  photoevents  per  image  is  small, 
and  the  covariance  of  X is  reasonably  constant  over  the  field  of  view, 
then 

^*0^  ' (3M 

The  final  problem  is  to  find  the  filter  matrix  [H]  that  minimizes 
the  error  measure  i of  Eq.(29)«  Fortunately,  the  answer  to  this 
problem  is  available  in  the  literature  and  can  be  directly  employed  here 
(Ref.l).  The  filter  matrix  that  minimizes  4 is  given  by 

•I 


(Hj  - k N(Sl(* J(B)‘| 


\ ^ 0 
• (Bjo 

_0  X 


♦ N(BH«^1(BJ‘)  . (3S) 


It  should  be  emphasized  that  this  filter  will  in  general  be  space  variant, 

depending  on  the  properties  of  (S]  , (B] . o , and  [t  ]. 

~ o 
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(c)  Minimum  Hean-Squared  Error 

In  order  to  specify  the  performance  of  the  optimum  filter,  It  is 

helpful  to  have  an  expression  for  the  minimum  mean-squared  error  4 t 

.mi  n 

actually  achieved.  If  we  substitute  expression  (35)  for  the  optimum 
fil'.er  in  expression  (2S)  for  the  mean-squared  error,  we  obtain 


^min 


Tr{tS][«j((l) 


- wH 


tB]o 


(36) 


where  (I]  is  the  identity  matrix.  Sometimes  it  is  useful  to  have 
the  minimum  mean-squared-error  in  the  alternative  form 

2 rr^ 


N Tr((S)(5ij(Sl^  - [Hj(B][ij(S)^)  (37) 


•min  \ 0 

For  purposes  of  comparison  with  the  continuous  case,  it  is 
interesting  to  manipulate  (36)  further,  yielding  the  form 


»min  • 


\ . 

0 

(BJo 

0 

X 

0 

[B)o 

0 

♦ N(B)  [*^HB)‘ 


-I 


[SlV 


(38) 


Note  the  similarity  of  this  result  to  expression  for  mln-mum  mean- 
squared  error  in  the  continuous  case  (derived  in  our  earlier  reports), 


(sec,  for  exemple,  E<i.(7)  of  RAOC-TR-77-165) . An  iinportent  difference 
Is,  of  course,  that  Eq.(38)  holds  In  the  space-variant  case,  while 
Eq.(39)  Is  valid  only  for  a space- Invariant  filter. 

(d)  Image  QMallty  Heasures 

Two  different  Image  quality  measures  were  found  useful  In  the 
continuous  case.  One,  which  we  called  the  "Image  quality  factor",  Q , 
was  defined  as  the  ratio  of  the  "signal"  energy  at  the  output  of  th.i 
restoration  filter  to  the  mean-squared  error  ("noise")  at  the  output. 

A second,  which  we  called  the  "restored  bandwidth",  tS  , was  equivalent 

e 

bandwidth  of  the  cascade  of  the  blur  and  deblurring  filters. 

Concentrating  first  on  the  quality  factor  Q , we  note  that  In  the 
discrete  case,  the  expected  energy  In  the  detected  signal  is  given  by 
the  Inner  product 


The  "signal"  portion  of  this  energy  is  given  by  the  second  term. 

After  passage  through  the  restoration  filter  (h]  , the  signal  energy  at 
the  output  Is  given  by 

Substitution  of  the  optimum  [h]  yields 
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ikl) 


k^ir  Tr|(S](*J(B]‘| 

r\  . 


[B)o 

.0  '\J 


-! 


(Bl  (*J[Bl 


(B]o 

0 \ 


♦ N[B][*JlBlV 


1-1 


♦N(B][*JfB]M  tBH*J[Sj' 


ik3) 


The  image  quality  factor  Q is  then  defined  as 

Q - m 

*min 

As  for  the  second  measure  of  image  quality,  restored  bandwidth, 
the  situation  is  a bit  more  difficult,  since  for  a space-variant  system 
there  is  no  unique  bandwidth  that  can  be  specified.  Stated  another  way, 
there  is  in  general  a different  bandwidth  associated  with  every  point 
in  the  f i I tered  image 

In  the  continuous  analysis  of  a space- invar iant  system,  the  restored 
bandwidth  aIT  was  defined  by 

- II  !BH|d»;jjdfly  (i»5) 

If  the  product  of  B and  M (the  normalized  transfer  functions  of  the 
blur  and  deblur  filters,  respectively)  is  non-negative  a.id  real,  as  is 
generally  the  case  in  practice,  then  £q.(^5)  can  be  rewritten 

ter 

. [[  h(-C.-n)bk.n)dtdn  (W) 

h 

where  b and  h are  the  inverse  Fourier  transforms  of  B and  H , 
respectively.  This  expression  is  the  integral  of  the  product  of  the 
two  impulse  responses,  centered  at  the  origin  in  this  case.  The  same 
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restored  bandwidth  would  be  obtained  by  integrating  the  product  of  the 
impulse  responses  centered  at  any  point  in  output  plane,  a consequence 
of  the  fact  that  both  the  blur  and  the  restoration  filter  are  space 
invar iant . 

In  the  space-variant  case,  the  product  of  b and  h will  yield 
functions  with  different  areas,  depending  on  where  this  product  is 
centered  in  the  output  plane  (x,y).  The  integral  becomes 

n(An'(x,y))^  - II  h(x,yU,n)b(C»n;x,y)dCdn  ?*7 ) 

• CM 

and  Is  clearly  a function  of  which  output  coordinates  (x,y)  are  chosen. 

The  discrete  analogy  to  Eq.(^7)  would  be  to  define  the  restored 
bandwidth  at  the  n^^  output  point  by 

. ((H)(Bl)„„  (48) 

where  for  metrix  [A]  , the  symbol  ((A))  1$  the  n**'  diegonal 

nn 

element  of  [A],  Since  this  quantity  depends  on  the  point  n , It  is 
not  of  itself  a useful  measure  of  image  quality*  However,  all  three  of 
the  following  definitions  would  be  useful: 


gTr((HHBl) 

Maximum  diagonal  element  of 

(H](Bl 

'A3b) 

Minimum  diagonal  element  of 

IHI(B). 

(Ii9c) 

These  definitions  may  be  regarded  as  yielding  average,  maximum  and 
minimum  restored  bandwidchs,  respect i vely . 

(e)  The  Stationary,  Shi  ft- Invar iant  Case 

We  now  consider  the  form  taken  by  the  various  results  above  in  the 
special  case  of  a shif t-invar iani  matrix  (h).  In  this  case,  o is  a 
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takes  on 


I ^ 

constant  vector  with  elements  — » while  the  matrix  K ] 

H 0 

toeplitz  form.  In  addition,  the  ideal  filter  matrix  [S]  and  the  blur 
matrix  [B]  take  toepi i taf  form,  and  the  matrix 
the  form  g H 1 • 

If  the  covariance  function  of  the  object  and  the  point-*spread 
functions  of  the  ideal  and  blur  filters  are  sufficiently  narrow,  the 
toeplitz  matrices  above  can  be  replaced  by  circulant^ approximations. 
With  this  substitution  the  post -compensat ion  filter  matrix  becomes 

[Hj  - l(8j*j  ' (SO) 

where  the  subscript  c indicates  a circulant  approximation. 

The  quality  factor  Q becomes  in  this  case 


Tr| 

[IH  JIB  J(i  ][B  J*(HJ*] 

W w C w j 

k^Trj 

'jdl  ♦ NlBjli  1(6  . 

l"  ' °c  ' 

'■) 

f ill)  IV'] 

The  equivalent  bandwidth  becomes 

-.(AJ5av>^  ■ i '52^ 

For  any  circulant  matrix  [C],  the  identity 

IC]  - (51IAJIJ)''  (5^) 

holds,  where  (9]  is  the  two-dimensional  discrete  Fourier  transform 
(OFT)  matrix,  and  [a^]  is  a diagonal  matrix  with  elements  that  are  the 
eigenvalues  of  (C).  Using  this  fact,  the  [k]  matrix  and  the  quality 
measures  can  be  rewritten  in  the  forms 

(Hj  - k MtJjlAjllAj,  JlAgl‘((I]  . NlAgllAj  1 (A  )*)'' (J)*' 

o o 

(5*t) 

*See  the  Appendix  for  , det,iled  definition  of  this  term. 


\ . O' 

[B]o 

!0  \ 


will  be  of 
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(55) 


Tr([y]tA„][AB][Aj,  ][Ab1’'([IJ  + NlAgl  tA,  1 [AbJ")' ' If] ’ 
k^(Tr([f)[As][Aj,  ]({1)  + NlAgllAj,  KAbI"  )' lAjJ^Ifl'') 

ii(Al!^y)^  • ^ Tr((f](A„][ABJ  (>]■')  (56) 

Note  that  the  eigenvalues  of  a clrculant  matrix  are  the  values  of  the 
OFT  of  the  first  column  of  the  clrculant  matrix. 

Because  of  the  diagonal  forms  of  (5^ , (55)  and  (56),  it  Is 
convenient  to  represent  the  n^^  component  of  the  OFT  of  the  restored 
Image  by  R(n)  , and  the  n^^  component  of  the  OFT  of  the  detected 
Image  by  0(n).  Then  we  have 

k N S(n)B''(n)i  (n) 

R(n)  5-2 0(n)  (57) 

I + N(B(n)ft^(n) 

and 

2.  |B(n)|‘'|S(n)|^*^'') 

N ). 2 

" (1  + M|B(n)|S„(n)) 

Q L (58) 

j*  |s(n)  I 4>^(n) 

n 1 ^ N|B(n) l^^^(h) 

A V 

where  ♦ Is  the  DM*  of  JL  Further, 
o ” 

n(A»^y)^  ■ R I H(n)B(n)  (59) 

n 

These  test  three  equations  are  exact  analogs  of  the  results  presented 
in  earlier  reports  for  the  continuous  case,  the  only  difference  being 
the  replacement  of  integrals  by  discrete  cums.  Thus  the  results  of  the 
more  general  space-var lent  analysis  do  Indeed  reduce  to  the  results 
obtained  in  earlier  analyses  of  the  space- Invar  lent  case. 
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V.  NUMERICAL  COMPUTATIONS  IN  PROGRESS 


At  this  time  computations  are  being  performed  on  four  simple  initial 
cases.  These  initial  cases  all  involve  one  dimensional  objects  only, 
and  each  has  at  most  one  shift-variant  or  non-stat lonary  factor.  One 
case  is  totally  stationary  and  sh i ft- invar iant . 

For  each  of  these  cases  the  image  quality  factor  Q of  (44)  and  the 
three  forms  of  the  restored  bandwidth  (49)  are  being  computed.  In  the 
three  non-stat ionary  or  shift-variant  cases  the  image  quality  measures 
are  being  computed  using  a series  of  shift- invar iant  restoration  filters 
[H] . In  the  totally  stationary,  shift-invariant  case  only  the  optimum 
shi f t- i nvar iant  [H]  are  being  used.  For  the  three  non-stat lonary  or 
shift-variant  cases  the  image  quality  measures  are  also  being  computed 
using  the  optimum  shift-variant  filter  of  (35)»  The  shift-variant  and 
shift- invar iant  filter  performances  will  then  be  compared. 

Of  our  initial  cases,  case  I is  totally  stationary  and  shift- 

A 

invariant.  For  this  case  o is  assumed  to  be 

£ " gi  (60) 

where  2.  ^ vector  of  all  ones.  In  general  the  following  relation 

holds 

_ t 

l5lol  - (♦qI^oo  (61) 

where 

_ t 

(♦qI  ■ E{(o-o)(or£)  } (62) 

In  the  stationary  case  the  covariance  matrix  (4^]  is  assumed  to  be 

(♦J  ■ o2[I)  (63) 

For  this  shift- Invar iant  case  the  blur  is  assumed  to  produce  a Gaussian 
point-spread  function.  In  each  of  the  following  three  cases,  only 
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one  of  the  above  assumptions  about  o , , and  [b]  wM)  be  changed. 

Thus  It  should  be  possible  to  determine  what  factors  are  most  improved 
by  the  shift-variant  filter  [h]. 

In  case  II  the  blur  [b]  Is  no  longer  assumed  to  be  shift- invar iant , 
but  Instead  Is  assumed  to  have  a position  dependent  width  which 
Increases  near  the  edges  of  the  image.  All  other  factors  are  assumed 
the  same  as  case  I . 

In  case  Ml  the  object  ensemble  is  assumed  to  have  a Gaussian  mean 
intensity  profile  centered  In  the  field  of  view.  All  other  factors 
are  again  assumed  the  same  as  case  I. 

In  the  last  case,  case  IV,  the  auto-covariance  function  [(^^]  is 
assumed  to  be  of  the  form 


m 


w • a width  parameter 


(65) 


0 - center  variance 

Thus  it  is  assumed  that  the  object  ensemble  has  a non-stat ionary 
variance  which  decreases  toward  the  edges  of  the  object.  As  in  the  above 
tM>  cases,  the  other  factors  are  assumed  the  same  as  case  I. 

Results  of  these  calculations  will  be  presented  in  our  next  technical 
report . 
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APPENDIX 


In  this  appendix  we  define  several  specialized  terms  used  in  the 
body  of  the  report. 

Lexicographic  ordering 

It  Is  common  practice  In  the  Image  processing  community  to  represent 
two-dimensional  arrays  of  samples  by  vectors.  If,  for  example,  an  Image 
Is  sampled  on  a two-dimensional  sampling  grid.  It  Is  perhaps  most 
naturally  represented  by  a matrix  of  sample  values: 


'll 

'12 

IN 

'21 

CM 

....  12,^ 

'ni 

'n2 

— 'nn 

(A.l) 


For  mathematical  purposes.  It  Is  often  convenient  to  reduce  this  matrix 
to  a vector  J_.  This  is  accomplished  by  scanning  the  matrix  row  by  row 
to  create  a one-dimensional  array  of  samples: 


11 


I 

I 


IN 

21 


I 


2N 


The  vector  Is  said  to  have  been  constructed  by  lexicographic 
ordering  of  the  matrix  [I]. 


(A. 2) 
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Toep] } tz  Matrix 


A matrix  [A]  is  said  to  be  a Toepi itz  matrix  if  it  has  the  property 

that,  whenever  j - k ■ m - n , the  elements  [A],,  and  [A]  are 

J K mn 

equal.  An  example  of  a Toepi itz  matrix  is  given  below: 


^ 7 5 

2 7 (A. 

3 2k 

6 9 2. 

The  covariance  matrix  of  a wide-sense  stationary  random  process  is 
always  a Toepi itz  matrix. 

Circulant  Matrix 

A matrix  [C]  is  called  circulant  if  each  of  its  rows  is  a right 
cyclic  shift  of  the  previous  row  by  one  element.  The  first  row  must  be 
a right  cyclic  shift  of  the  last  row  by  one  clement.  An  example  of  a 
circulant  matrix  is  given  by 


IC]  - 


2 k 
5 2 
7 5 

^ 1 


7 5 

k 7 
2 k 
5 2 


(A.i») 


Note  that  by  "cyclic  shift"  we  mean  the  shifting  of  rows  must  be  of  the 
"wrap-around"  type.  That  is,  when  an  element  Is  shifted  out  on  the  right, 
it  immediately  reappears  on  the  left.  Note,  also  that  all  circulant 
matrices  are  Toepi itz  In  form. 

It  is  common  practice  in  image  processing  to  approximate  some  forms 
of  Toepi itz  matrices  by  constructing  circulant  matrices  from  them.  If 
a Toepi itz  matrix  has  only  a few  localized  non-zero  diagonals,  a circulant 
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approximation  Is  often  constructed  by  replacing  some  of  the  zero  diagonals 


with  non*zero  diagonals  thus  forming  the 
As  an  example,  the  Toeplitz  matrix 


[T] 


7 5 I 
3 7 5 
2 3 7 
0 2 3 
0 0 2 
0 0 0 
0 0 0 


can  be  converted  to  the  circulant  matrix 


IT'I 


7 5 I 
3 7 5 
2 3 7 
0 2 3 
0 0 2 
I 0 0 
5 I 0 


result  into  the  circulant  form. 


0 0 0 0 

10  0 0 
3 10  0 

7 5 10 

3 7 5 1 

2 3 7 5 

0 2 3 7 


(A. 5) 


0 0 2 3 
10  0 2 
5 10  0 
7 5 10 
3 7 5 1 
2 3 7 5 
0 2 3 7 


(A.6) 


by  changing  the  three  upper  right  and  lower  left  zero  elements. 
Computations  using  circulant  matrices  can  be  performed  faster  than  those 
for  Toeplitz  matrices  because  the  discrete  Fourier  transform  can  be 
utilized.  As  the  fraction  of  diagonals  of  the  Toeplitz  matrix  which  are 
non-zero  decreases,  the  circulant  approximation  approaches  the  original 
matrix,  and  the  errors  associated  with  the  use  of  the  circulant  approxi- 
mation can  be  shown  to  decrease. 
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MISSION 

of 

Rome  Air  Development  Center 


RADC  plans  and  conducts  research,  exploratory  and  a&iranced 
development  programs  in  coamand,  control,  and  coammications 
(C^)  activities , and  in  the  areas  of  informatioTi  sciences 
and  intelligence.  The  principal  technical  mission  areas 
are  conmunications , electromagnetic  guidance  and  control, 
surveillance  of  ground  and  aerospace  objects,  intelligence 
data  collection  imd  handling,  information  system  technology, 
ionospheric  propagation,  solid  state  sciences,  siicrotaave 
physics  and  electronic  reliability , maintainability  and 
compatibility . 


